%% Simulation with minimax
% probability alpha = 0->1
% C = E[alpha]

clc;
clear all;
close all;
%% Initialization
N = 4;
m = N;n=N;r=N;
k = 1;

alpha = 0:0.01:1;%prior

%step of multiplicand, multiplier and correction value
step_a = 2^(-m);
step_b = 2^(-n);
step_c = 2^(-r-k);
ulp = 2^(-r);

%range of a,b
lower = 0.0;
upper = 0.999999;

num_a = 2^m;
num_b = 2^n;
num_ab = num_a * num_b;

mae = zeros(1, length(alpha));%maximum absolute error
mse = zeros(1, length(alpha));%expected absolute error
ave = zeros(1, length(alpha));%average error
for i = 1:length(alpha)
    pi = alpha(i);%prior
    cvi = ccm_rd(m,n,r,k,pi);%expected value of reduction error
    cvi_r = rnd_near(cvi,r+k);
    for a = lower:step_a:upper
        for b = lower:step_b:upper
            z = a * b;%true product
            z_r = rnd_near(z,r);%n bit product
            
            z_pi = trunc_mult(a,b,cvi_r,m,n,r,k);
            err = z_pi - z;
            abs_err = abs(err);
            if(mae(i) < abs_err)
                mae(i) = abs_err;
            end
            
            mse(i) = mse(i) + abs_err^2;
            ave(i) = ave(i) + err;                        
        end
    end
    mse(i) = mse(i)./num_ab;
    ave(i) = ave(i)./num_ab;
end

%% Plot
close all;
plot(alpha, mae./ulp);
figure;
plot(alpha, mse);
figure;
plot(alpha, ave);
